#setup: install packages if not installed already
installPackageNotFound <- function(package_name){

  possibleError <- tryCatch(
    suppressPackageStartupMessages(library(package_name, character.only = TRUE)),
    error=function(e) e
  )
  if(inherits(possibleError, "error")){
    install.packages(package_name, repos=c("http://cran.stat.ucla.edu",
    "http://cran.rstudio.com/"))
    suppressPackageStartupMessages(library(package_name, character.only = TRUE))
  }

}

## Function to format number to a specific number of sig digits
formatSig <- function(number, n){
  return(formatC(round(number, n), digits = n, format="f", flag="#"))
}

#Function to run t-tests split by a population variable for a particular outcome
#and independent variable
#Applies a Benjamini Hochberg correction to p values within population
#splits to adjust for multiple tests
runTTests <- function(test.data, outcome.var.names, indp.var){
  #############################
  #Keyword args
  #test.data = dataframe containing columns named for outcome.var and indp.var
  #outcome.var.names = character vector of variable names to use for outcome
  #indp.var = string name of variable to use for independent variable

  #Returns
  #dataframe with the test results
  #############################

  #Prep loop
  results.df <- c()

  #Run t-tests for each covariate
    for(outcome.var in outcome.var.names){
      #try t-test
      possibleError <- tryCatch(
        ttest <- t.test(test.data[[outcome.var]]~test.data[[indp.var]]),
        error=function(e) e
      )
      if(inherits(possibleError, "error")){
        tcount <- table(test.data[[indp.var]])
        test.row <- c(outcome.var,NA,NA,NA,NA,tcount[1],tcount[2], NA, NA)
      }else{
        sdev <- aggregate(formula(paste0(outcome.var,"~",indp.var)),
        data = test.data, FUN = sd)
        tcount <- table(test.data[[indp.var]])
        test.row <- c(outcome.var,ttest$estimate[1], ttest$estimate[2],
                      sdev[1,outcome.var]/sqrt(tcount[1]),
                      sdev[2,outcome.var]/sqrt(tcount[2]),
                      ttest$statistic,ttest$p.value, tcount[1], tcount[2],
                      ttest$parameter[['df']],
                      (ttest$estimate[1]-ttest$estimate[2])/ttest$statistic)
      }
      results.df <- rbind(results.df, test.row)
    }

  #Clean up results
  results.df <- suppressWarnings(data.frame(results.df, stringsAsFactors = F))
  names(results.df) <- c("outcome", "control_mean", "treatment_mean",
                         "control_se", "treatment_se",
                         "t_value", "p_value", "n_control", "n_treatment",
                         "df", "se.diff")

  #Convert numbers to numeric
  results.df[,c(2:ncol(results.df))] <- sapply(results.df[,c(2:ncol(results.df))],
  as.numeric)

  #Adjust p-value
  results.df$bh_p_value <- p.adjust(results.df$p_value, method="fdr")

  #Add column for difference
  results.df$difference <- results.df$treatment_mean - results.df$control_mean

  #Add stars
  results.df$sig.stars <- ""
  results.df$sig.stars[which(results.df$bh_p_value < 0.001)] <- "***"
  results.df$sig.stars[which(results.df$bh_p_value >= 0.001 &
    results.df$bh_p_value < 0.01)] <- "**"
  results.df$sig.stars[which(results.df$bh_p_value >= 0.01 &
    results.df$bh_p_value < 0.05)] <- "*"

  return(results.df)
}
